%%%% Solution Algorithm for the Representative-Agent (with Endogenous Labor) Equivalent for the Aiyagari (1994,QJE)                  
%%%% Written by Orhan Torul & Oguz Oztunali
%%%% October 2017, Bogazici University, Istanbul

%%
%%
%% 0. housekeeping
clc;                                % clear screen
clear;                              % clear memory
close all;                          % close open windows
%%
%% 1. parameters and functional forms

% parameters
alpha  = 0.56;                      % capital income share from Penn World Table 9.0 (PWT)
b      = 0;                         % borrowing constraint as in Aiyagari (1994, QJE)
beta   = 89/100;                    % subjective discount factor from PWT
delta  = 5.5/100;                   % depreciation rate of physical capital from PWT
gamma  = 3/2;                       % inverse elasticity of intertemporal substitution from macroeconomics literature
varphi = 2/3;                       % Frisch elasticity of labor supply from macroeconomics literature
zbar   = 0.783472854785600;         % effective labor productivity from the heterogeneous-agent incomplete market model


% utility and (inverse) marginal utility functions
% contemporaneos utility is defined as u(c)-v(x) with the functinal forms below

u     = @(c) (c.^(1-gamma))/(1-gamma);          % utility from consumption
v     = @(h) (h.^(1+1/varphi))/(1+1/varphi);    % disutility from labor supply
up    = @(c) c.^(-gamma);                       % marginal utility of consumption
invup = @(x) x.^(-1/gamma);                     % inverse of marginal utility of consumption
vp    = @(h) h.^(1/varphi);                     % marginal disutility of labor supply
invvp = @(x) x.^(varphi);                       % inverse marginal disutility of labor supply 



%% 2. deterministic steady-state
 
% initial guess
y0 = [10,0.5,1,2,0.12]; 
% assignment of endogenous variables are as follows:
% y(1) capital (k)
% y(2) labor (h)
% y(3) consumption (c)
% y(4) wage (w)
% y(5) real interest rate (r)

% option setting for the fsolve function
options= optimset('maxiter',999999999999','tolfun',10e-8,'tolx',10e-8);
[SS, fval,exitflag,output]=fsolve(@(y) steady(y,alpha,beta,delta,gamma,varphi,zbar),y0,options);

kss=SS(1);
hss=SS(2);
lss=hss*zbar;
css=SS(3);
wss=SS(4);
rss=SS(5);
gdpss=kss^alpha*lss^(1-alpha);
incomess=gdpss-delta*kss;

%% 3. reporting results

disp ('steady-state for the decentralized economy:')
disp('-------------------------------------------------')
fprintf('output (y):              %.4f \n',gdpss);
fprintf('capital (k):             %.4f \n',kss);
fprintf('labor (h):               %.4f \n',hss);
fprintf('effective labor (l):     %.4f \n',lss);
fprintf('consumption (c):         %.4f \n',css);
fprintf('wage (w):                %.4f \n',wss);
fprintf('real interest rate (r):  %.4f \n',rss);
disp('-------------------------------------------------')
fprintf('utility from c:          %.4f \n',u(css));
fprintf('disutility from n:       %.4f \n',v(hss));
fprintf('contemporaneous utility: %.4f \n',u(css)-v(hss));

meanmatrix=[kss hss zbar lss gdpss css rss wss rss*kss wss*lss ];

% exporting steady-state values to LaTeX
Det_TableStationaryMeans=[meanmatrix];
columnlabelStationaryMeans = {'$K$', '$H$' ,'$Z$' ,'$L$', '$Y$', '$C$', '$r$', '$w$' ,'$rK$', '$wL$'};
matrix2latex(Det_TableStationaryMeans, 'Det_TableStationaryMeans.tex', 'columnLabels', columnlabelStationaryMeans,'format', '%-6.3f', 'size', 'footnotesize');